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Abstract 

We report a theoretical study of a quantum optical model consisting of an array of strongly nonlinear cavities inco¬ 
herently pumped by an ensemble of population-inverted two-level atoms. Projective methods are used to eliminate 
the atomic dynamics and write a generalized master equation for the photonic degrees of freedom only, where the 
frequency-dependence of gain introduces non-Markovian features. In the simplest single cavity configuration, this 
pumping scheme gives novel optical bistability effects and allows for the selective generation of Fock states with a 
well-defined photon number. For many cavities in a weakly non-Markovian limit, the non-equilibrium steady state 
recovers a Grand-Canonical statistical ensemble at a temperature determined by the effective atomic linewidth. For a 
two-cavity system in the strongly nonlinear regime, signatures of a Mott state with one photon per cavity are found. 
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1. Introduction 


The study of quantum many-body systems is one of the most active fields of modern condensed-matter physics. 
Among the most celebrated effects, we can mention frictionless flows in superfluid and superconducting systems 
and the geometrical quantization features of the fractional quantum Hall effect. While this physics was traditionally 



in atomic nuclei @], in quark-gluon plasmas AS, or in electron gases 
the last two decades have witnessed impressive advances using ultra-cold 


studied in liquid Helium samples 
confined in solid-state devices [6,(7 

atomic gases trapped in magnetic or optical traps JT3UHL d. 

In the last few years, a growing community has started investigating many-body effects in the novel context of 
the so-called quantum fluids of light 11311 . i.e. assemblies of many photons confined in suitable optical devices, 
where effective photon-photon interactions arise from the optical nonlinearity of the medium. After the pioneering 
studies of Bose-Einstein condensation Cl and superfluidity Cl effects in dilute photon gases in weakly nonlinear 
media, a great interest is presently being devoted to strongly nonlinear systems, where even single photons are able to 
appreciably affect the optical properties of the system. 

The most celebrated example of such physics is the photon blockade effect d, where the presence of a single 
photon in a cavity is able to detune the cavity frequency away from the pump laser, so that photons behave as ef¬ 
fectively impenetrable particles. Experimental realizations of this idea have been reported by several groups using 
very different material platforms, from single atoms in macroscopic cavities d, to single quantum dots in photonic 
crystal cavities 11_8, 19j], to single Josephson qubits in circuit QED devices for microwaves |2Cj, 21]. 

Scaling up to arrays of many cavities coupled by photon tunneling is presently a hot challenge in experimental 
physics, as it would realize a Bose-Hubbard model for photons where the photon blockade effect may lead to a rich 
physics, including the superfluid to Mott-insulator phase transition at a commensurate filling or Tonks-Girardeau 
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gases of impenetrable photons in one-dimensional continuum models. The first works on strongly correlated photons 
were restricted to quasi-equilibrium regimes where the photon loss rate is much slower than the internal dynamics of 
the gas so that the system has time to thermalize and/or be adiabatically transfered to the desired strongly correlated 
state [22j,l23]. While this assumption might be satisfied in suitably designed circuit-QED devices in the microwave 
domain, radiative losses are hardly negligible in realistic optical cavities in the infrared or visible domain, so that 
thermalization is generally far from being granted | T31 2l|]. 

As a result, a very active attention has been recently devoted to the peculiar non-equilibrium effects that arise for 
realistic loss rates. Starting from the pioneering work on photon blockade in non-equilibrium photonic Josephson 
junctions {241]. the interest has been focused on the study of schemes to generate strongly correlated many-body 
states in the very non-equilibrium context of photon systems, where the steady-state is not determined by a thermal 
equilibrium condition, but by a dynamical balance of driving and losses. 

The first such scheme proposed in [23] was based on a coherent pumping: provided the different many-body 
states are sufficiently separated in energy, many-photon processes driven by the coherent external laser are able to 
selectively address each many-body state as done in optical spectroscopy of atomic levels. In this way, the non¬ 
equilibrium condition is no longer just a hindrance, but offers new perspectives, as it allows to individually probe 
each excited state. Furthermore, the appreciable radiative losses make microscopic information on the many-body 
wavefunction be directly encoded in the quantum coherence of the secondary emission from the device 2fJ 2J, 281. 
While this coherent pumping scheme offers a viable way to generate and control few photon states in small arrays, 
its efficiency is restricted to mesoscopic systems where the different states are well-separated in energy. Moreover, 
this scheme intrinsically leads to coherent superpositions of states of different photon number: while this feature 
is intriguing in view of observing many-body braiding phases {2811 , it is not ideally suited to generate states with a 
well-defined photon number such as Mott-insulator states. 

The identification of new schemes that do not suffer from these limitations is therefore of great importance in 
view of experiments. In the present work we study the potential of frequency-dependent gain processes to selectively 
generate strongly correlated states of photons in arrays of strongly nonlinear cavities. The frequency-dependence 
of amplification is a well-known fact of laser physics and is often exploited to choose and stabilize a desired lasing 
mode [29]. In the last years, a series of works by our groups I 30l 31 1 have explored its effect on exciton-polariton 
Bose-Einstein condensation experiments, in particular questioning the apparent thermalization of the non-condensed 
fraction ^32, 33, 313 S’]- All these works were however restricted to the weakly interacting regime where quantum 
fluctuations can be treated in the input-output language by means of a Bogoliubov-like linearized theory around 
the mean-field. Here we tackle the far more difficult case of strong nonlinearities, which requires including the 
non-Markovian features due to the frequency-dependent gain into the many-body master equation for the strongly 
interacting photons and then to solve the quantum many-body theory of the generalized driven-dissipative Bose- 
Hubbard model. 

In the last years, similar questions have been theoretically addressed by several groups. Just to mention a few of 
them, a scheme to obtain a thermal state at finite temperature with a non-vanishing effective chemical potential for 
photons has been proposed in [ 3711 using a clever parametric system-bath coupling with as special eye to circuit-QED 
and opto-mechanical systems. A further development in this direction 113811 has considered pumping by two-photon 
processes in the presence of an auxiliary shadow lattice in a circuit-QED architecture: in spite of the complexity of 
the proposed set-up, the mechanism underlying the stabilization of many-body states is very similar to our frequency- 
dependent gain. With respect to these proposals and to the engineered dissipations originally proposed for atoms 13911 
and then extended to photons [40. 4l| to organic polaritons 142. 43;], and circuit QED systems [44, 3 5J, our 
approach has the crucial advantage of being based on a quite commonly observed feature of laser and photonic 
systems such as a frequency-dependent gain. Finally, a pioneering discussion of the onset of collective coherence in a 
related model of a cavity array embedding population-inverted atoms has recently appeared in [471], but little attention 
was paid to the effect of strong nonlinearities nor to the development of a tractable quantum formalism. 

The aim of this article is to introduce the readers to the basic physics of a frequency-dependent incoherent pump¬ 
ing and to first illustrate the consequences of the resulting non-Markovianity in the simplest configurations before 
attacking more complex many-body effects. With this idea in mind, the structure of the article is the following. In 
Sec [2] we present the physical system and we develop the theoretical model based on a master equation for the cavities 
coupled to the atoms of the gain medium. The projective method to eliminate the atomic degrees of freedom and write 
a master equation for the photonic density matrix is sketched in Sec l2.2l along the lines of the general theory of [|48]. 
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First application of the method to a single cavity configuration is discussed in Sec|3]and specific features of the weak 
and the strong nonlinearity cases are illustrated, e.g. a novel mechanism for optical bistability and the selective gener¬ 
ation of Fock states with a well defined photon number. The richer physics of many cavity arrays is discussed in Sec(4j 
In a Markovian regime, the photonic steady state has the surprisingly trivial form of a Grand-Canonical distribution 
of infinite temperature, and therefore is fully independent of the many-body photonic Hamiltonian. In a weakly non- 
Markovian regime, an effective Grand-Canonical distribution of finite temperature is obtained even in the absence of 
thermalization mechanisms; in a strongly nonlinear and non-Markovian regime, signatures of a Mott insulator state 
with one photon per cavity are illustrated. Conclusions are finally drawn in Sec| 6 ] In the Appendices, we provide 
the details of the derivation of the photonic master equation using projective methods, on the exact stationary state in 
the Markovian case, on a perturbative expansion of the coherences in the weakly non-Markovian limit, and on further 
numerical validation of the purely photonic master equation. 


2. The physical system and the theoretical model 


2.1. The physical system 

In this work, we consider a driven-dissipative Bose-Hubbard model for photons in an array of k coupled nonlinear 
cavities of natural frequency uj cav . In units such that h = 1, the Hamiltonian for the isolated system dynamics has the 
usual form 113, 21 .j49ll: 


k 

Hp h = 

t ^ f f 
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(1) 


They are arranged in a one-dimensional geometry and are coupled via tunneling processes with amplitude J. Each 
cavity is assumed to contain a Kerr nonlinear medium, which induces effective repulsive interactions between photons 
in the same cavity with an interaction constant U proportional to the Kerr nonlinearity . Dissipative phenomena 
due the finite transparency of the mirrors and absorption by the cavity material are responsible for a finite lifetime of 
photons, which naturally decay at a rate F; oss . 

As mentioned in the introduction, the key novelty of this work with respect to earlier work consists in the different 
mechanism that is proposed to compensate for losses and replenish the photon population. Instead of a coherent 
pumping or a very broad-band amplifying laser medium, we consider a configuration where a set of N at two level 
atoms is present in each cavity. Each atom is strongly pumped at a rate r pump , spontaneously decays to its ground 
state at a rate 7 and, most importantly, is coupled to the cavity with a Rabi frequency f2#: as a result, the atoms 
provide an incoherent pumping of the cavities, with a frequency-dependent rate centered at the atomic frequency ui a t- 
Our choice of two different physical mechanisms for nonlinearity and pumping (for example, two different atomic 
species) allows us to to tune independently photonic interactions and emission. 

The free evolution of the atoms and their coupling to the cavities are described by the following Hamiltonian 
terms, 

k N at 

ff„. = £5>«7 < ‘A < ‘ > 

2=1 Z = 1 

k N at 

Hl = ^EE[ a k ' (0 +^ +(0 
2—1 1=1 

the atomic frequency uj at is assumed to be in the vicinity (but not necessarily resonant) with the cavity mode and 
the atom-cavity coupling is assumed to be weak enough Un <C uiat, to cav to be far from the ultra-strong coupling 
regime ll50h and from any superradiant Dicke transition EH. 

As usual, the dissipative dynamics under the effect of the pumping and decay processes can be described in terms 
of a master equation for the density matrix p of the whole atom-cavity system, 

dtp = — [Hph + H at + Hi, p] + C(p), 
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where the different dissipative processes are summarized in the Lindblad super-operator C = C pump + C/ oss . at + 
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describing the pumping of the atoms, the spontaneous decay of the atoms, and the photon losses, respectively. The 
erf 1 ' 1 ' 1 operators are the usual raising and lowering operators for the /-th atom in the /-th cavity. We introduce the 
detuning 5 = u> cav — w a t of the bare cavity frequency with respect to the atomic frequency. In the following, we shall 
concentrate on a regime in which pumping of the atoms is much faster than their spontaneous decay, r pump 7 , so 
the Cioss, at Lindblad term can be safely neglected. 


For simplicity, we will also restrict our attention to the F 


pump 


• regime, where the atoms are immedi¬ 
ately repumped to their excited state after emitting a photon into the cavity: under such an assumption, an atom having 
decayed to the ground state does not have the time to reabsorb any photon before being repumped to its excited state. 
In this regime, complex cavity-QED effects such as Rabi oscillations do not take place and the photon emission takes 
place in an effectively irreversible way |48, 52j|: as a result, we are allowed to eliminate the atomic dynamics from the 


problem and write a much simpler photonic master equation involving only the cavity degrees of freedom. 


2.2. Closed master equation for the photonic density matrix 

Under the considered r pump 1 1 /,• approximation, the atomic population is concentrated in the excited state and 
it is possible to use projective methods to write a closed master equation for the photonic density matrix where the 
atomic degrees of freedom B have been traced out, p p h = Trgp. All details of the (quite cumbersome) calculations 
can be found in Appendix A The resulting photonic master equation reads 

tdtPph = t [fTp/i, Pphij)] T Cl oss T ^em 1 (8) 

with 
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describing photonic losses and emission processes, respectively. While the loss term has a standard Lindblad form at 
rate l’/ os .s, the emission term keeps some memory of the atomic dynamics as it involves modified lowering and raising 
operators 


a = a* 


n OO 

■ / dr e ( '~ iu ’ at ~ Tpump/2)T di(-T), 
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which contain the photonic (hamiltonian and dissipative) dynamics during pumping. In the limit we are considering 
in which photonic losses are slow with respect to atomic pump, these operators are the interaction picture ones with 
respect to the photonic hamiltonian in the cavity array and have a simpler expression : 

Oi(r) =e iH * hT a i e- iH ‘’ hT . (13) 
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The Fourier-like integral in Eqs[TT]and[T2]is responsible for the frequency selectivity of the emission, as the integral 
is maximum when the free evolution of a, occurs at a frequency close to the atomic one ui a t- 

A deeper physical insight on the operators ( ITTb and (IT2l) can be obtained by looking at their matrix elements 
in the basis of eigenstates of the photonic hamiltonian. We consider two eigenstates |/) (resp. |/')) with N (resp. 
N + 1) photons and energy w/ (resp. uip). After elementary manipulation, we see that the emission amplitude 
follows a Lorentzian law as a function of the detuning between the frequency difference of the two photonic states 
uif/f = ujfi — ojf and the atomic transition frequency ui at . 
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Upon insertion of Eq[14]into the master equation Eq|8j one can associate the real part of the Lorentzian factor to an 
effective emission rate 


r em(w/'/) — r em 


r 2 

pump 


/4 


iyJat ^pump 
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(15) 


while the imaginary part can be related to a frequency shift of the photonic states under the effect of the population- 
inverted atoms. In the next section, this point will be made more precise under a secular approximation. 

The width of the Lorentzian is set by the pumping rate r pump , that is by the autocorrelation time T pvrnp = 
1 /r pump of the atom seen as a frequency-dependent emission bath. The peak emission rate exactly on resonance is 
equal to 


T 


0 

em 


4,N at n 2 R 


- pump 


(16) 


While the T pump y/N a t£ln assumption automatically implies that the emission is much slower than the atomic 
repumping rate, Y em <C T pump , no constraint need being imposed on the parameters J, U and S = ui cav — uj a t of 
the photonic Hamiltonian, which can be arbitrarily large. Whereas an extension of our study to the r; oss > T pump 
regime would only introduce technical complications, entering the T em > T purnp regime is expected to dramatically 
modify the physics, as a single atom could exchange photons with the cavity at such a fast rate that it has not time 
to be repumped to the excited state in between two emission events. As a result, reabsorption processes and Rabi 
oscillations are possible, which considerably complicate the theoretical description. These issues will be the subject 
of future investigations. 


2.3. Reformulation in Lindbladform in the secular approximation 

In the case the system has a discrete spectrum, it is possible in the so-called secular approximation to write another 
photonic master implementing non-markovian effects with a more standard Lindblad form, compatible with Monte 
Carlo wave-function simulations [531] and giving equivalent driven-dissipative dynamics. This can be explained by 
the following argument: in a weak dissipation limit (T em , Tj oss very small with respect to the gaps in the spectrum) 
terms of the density matrix pj p pj, p which would be rotating at different frequencies uif p j, if the system 
were isolated, are not coupled to each other by dissipation since the coupling T° em , Ti oss is negligible with respect to 
their frequency difference Acu = oj f/ p — u> f f = topj — to f, f. Considering this, all relevant dissipative transitions 
verify then Aw ~ 0. Restricting the previous master equation given by Eas. 181 1 Ol and ll~4l to these transitions, it is 
possible to rewrite the dynamics in the following way (details of the derivation are given in | Appendix B 1 : 
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Note that the jump operators a] have the same form as the ones considered in [38] and have for effect to modify the 
the transition rate, while the "imaginary part" of Eq.[14]induces an additional Hamiltonian contribution in the form of 
a Lamb shift. Notice that the two master equations Eqs. (8j[T7]are slightly different. However, under the considered 
approximation they are expected to provide equivalent dynamics. The latter form has the advantage of being of 
Lindblad form, and thus is directly compatible with MCWF simulations [53] and can be useful from a numerical 
point of view. 

The secular approximation can be very restrictive (particularly in the thermodynamic limit where the spectrum is 
continuous). However, our feeling is that the reformulation of Eq|T7]should be accurate in a wider range of parameters. 
Quantitatively, we anticipate the condition T em , Ti oss -C T pum p to be sufficient. More investigations in this direction 
are under way. 


3. One cavity 

As a first example of application, we consider the simplest case of a single nonlinear cavity. A special attention 
will be paid to the stationary state p ss of the system for which Eq|8]imposes 


0 — 2 \Hph ) Pss] T t^loss (P ss ) + C em (/^ss ) • (21) 

In our specific case of a single cavity, the photonic states are labelled by the photon number N and have an energy 

UJN = Nco cav + ^N(N - l)U. (22) 

Correspondingly, the N —> N + 1 transition has a frequency 


d>N+l,N — Wcav + NU, 


(23) 


and the corresponding photon emission rate is 


r e m(wAr+l,Jv) = r 
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em 


{Tpump/^f 


(tUjV+l,tV T {Fpump/ 2)^ 


(24) 


As no coherence can exist between states with different photon number N, the stationary density matrix is diagonal 
in the Fock basis, p ss = Sn,n'^n with the populations 7r/v satisfying 


(N + l)r toss 7rjv + i — (N + l)r em (u;jv+i,iv)7rjv + NT ern (LON,N-i)^N-i — NTi oss ttn = 0, (25) 


where the two last terms of course vanish for N = 0. As only states with neighboring N are connected by the 
emission/loss processes, detailed balance is automatically enforced in the stationary state, which imposes the simple 
condition on the populations, 


(N + l)r; oss 7Tjv + i — (N + l)T em (oj N+1 , N )nN = 0 

which is straightforwardly solved in terms of a product, 


T/V = 7h) 
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The excellent agreement of this result with a numerical solution of the full atom-cavity system is illustrated in 


Appendix E 
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Figure 1: (a) Emission vs. loss rate as a function of the detuning from the atomic frequency u) a t : the three curves are for peak emission Tg m 
larger (red dash-dotted), equal (black dashed), smaller (green solid) than the loss rate Ti oss . (b-d) Populations ttn of the ./V-photon state as a 
function of N in the three cases LJ 2 < cu C av (b), oji < co C av < ^2 (c), cj ca ,; < cji (d). In the three panels, the open dots are the numerical 
results of the atom-cavity theory, while the solid line is the prediction of the analytical purely photonic theory; the dashed curves show the ratio 
Fem(wN+i,N)/Floss as a function of N. Parameters: 8/U = 4 (b), —2 (c), —6 (d). In all panels, 2U/T pU mp = 0.2, 2Ti oss /V pU mp = 
0.0006, 2£ln/Y purnp = 0.02. 
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3.1. Linear regime 

For a vanishing nonlinearity U = 0, all transition frequencies un+i,n are equal to the bare cavity frequency uiq 
and the populations of the different N states have a constant ratio 


7r A f +1 _ r( Lpump /2) 

if n r; oss (5 2 + (T pump /2) 2 ’ 

where we remind that S = oj cav — ui a t- For weak pumping and/or large detuning, one has 


(Tpump/2) 2 

l S 2 + {Tpump/2/ 


<r 


loss f 
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so the density matrix for the cavity shows a monotonically decreasing thermal occupation law. 
and close to resonance, one can achieve the regime where the emission overcompensates losses 
starts being strongly populated: 


(Tpump/2) 2 

'S 2 + (T pU m P /2) 


>r 


loss • 


For strong pumping 
and the cavity mode 

(30) 


The transition between the two regimes is the usual laser threshold, but our purely photonic theory is not able to 
include the gain saturation mechanism that serves to stabilize laser oscillation above threhsold | 29| , 521: within our 
purely photonic theory, the population would in fact show a clearly unphysical monotonic growth for increasing N. A 
complete description in terms of the full atom-cavity master equation would of course solve this pathology including 
a gain saturation mechanism according to usual laser theory, but this goes beyond the scope of the present work. 


3.2. Optical bistability phenomena in weak nonlinear cavities 

For U > 0, the situation is much more interesting as the effective transition frequency depends on the number of 
photons, 

t^jV+ljjV — ktcav T NU ^ id cav , , (31) 

so the gain condition 

Tern _ (Tpump/2) 2 _ ^ ^ ^2) 

T/oss (*tN+l,N tU at ) 2 T (Tpump/2) 2 

can be satisfied in a finite range of photon numbers only, as it is illustrated in FigUfa). As a consequence, even a 
weak nonlinearity U is able to stabilize the system for any value of \' { / nj even in the absence of any gain saturation 
mechanism. 

For F° m < F; oss , losses always dominate. For T/, m > r; oss , the gain condition is instead satisfied in a range 
of frequencies [u>i, W 2 ] around io at . Under the weak nonlinearity condition U <C T pU mp , the [uji , W 2 ] range typically 
contains a large number of transition frequencies u>n+i,n at different N. Three different regimes can then be identified 
depending on the position of the cavity frequency uj cav with respect to the [wi, UI 2 ) range. 

(i) If w 2 < uj cav , then the gain condition is never verified, and the population tty shown in FigJTJb) is a mono¬ 
tonically decreasing function of N. In this regime, the state of the cavity field is very similar to a thermal state, as 
it usually happens in a laser below threshold, (ii) If w 1 < uj cav < 0 J 2 , the population shown in FigJTfc) is an 
increasing function for small N, shows a single maximum for N ~ N = ( 0 J 2 — u} cav )/U, and finally monotonically 
decreases for N > N. 

The phenomenology is the richest in the regime (iii) where uj cav < u>i. In this case, for small N the population 
ttn decreases from its initial value 7r (J until the nonlinearly shifted frequency enters in the gain interval for N ~ 
N’ = (w 1 — ui cav )/U. After this point ttn starts increasing again until it reaches a local maximum at N ~ N = 
(uj 2 — uJcav) /U. Finally, for even larger N it begins to monotonically decrease. An example of this complicate 
behaviour is shown in FigdJd). 

The existence of two well separate local maxima at N = 0 and N ~ N in the photon number distribution 
7 t/v suggests that the incoherently driven nonlinear cavity exhibits a sort of bistable behaviour: when it is prepared 
at one maximum of the photon number distribution 7 r/y, the system is trapped in a metastable state localized in a 
neighborhood of this maximum for a macroscopically long time. Switching from one metastable state to the other 
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Figure 2: Purely photonic simulation of the two-time coherence function g (r) in the weakly nonlinear regime. Parameters U/T p 
Fioss/^pump = 0.03, r2 m /r pump = 0.04, S = -6 U as in Fig[Qd). 


= 0.1, 


results is only possible as a result of a large fluctuation, so it has a very low probability, typically exponentially small 
in the photon number difference between the two metastable states. 

This bistable behavior is clearly visible in the temporal dependence of the delayed two-photon correlation function 


. 9 ( 2 ) 0 ) = 


(<zt(t) at(t + t) a{t + r) a(t)) ss 
(at(f) a(t)) ss (at(f + r) a(t + t)) s 


(33) 


that is plotted in Figj2] At short times, the value of g <2> is determined by a weighted average of the contribution of 
the two maxima according to the stationary ttn- After a quick transient of order l/r em / oss , which corresponds to a 
fast local equilibration of the probability distribution around each of its maxima, the g l2> correlation function slowly 
decays to its asymptotic value 1 on a much longer time-scale mainly set by the exponentially long switching time 
from one maximum to the other 

Before proceeding, it is worth emphasizing that the present mechanism for optical bistability bears important 
differences from the dispersive or absorptive optical bistability phenomena discussed in textbooks 1154115511. On one 
hand there is some analogy to dispersive optical bistability in that the intensity-dependence of the refractive index is 
responsible for a frequency shift of the cavity resonance; on the other hand the frequency-selection is not provided by 
the resonance condition with a monochromatic coherent incident field rather by the frequency dependence of the gain 
due to the incoherent pump. 


3.3. Photon number selection in strongly nonlinear cavities 

In the opposite limit U T pump , the nonlinearity is so large that a change of photon number by a single unity 
has a sizable effect on the emission rate r em (wjv+i,jv)- As discussed in Appendix A| the derivation of the photonic 
master equation remains fully valid in this regime provided T pum p ^ Tg m , r; oss . 

The ensuing physics is most clear in the regime when the maximum emission rate is large but only a single 
transition fits within the emission lineshape: these assumptions are equivalent to imposing that 


T° 


» 1 


and 


r° r 2 

L em pump ^ 


Tioss Ti oss U 2 

with the further condition that the emission is resonant with the Nq —> Nq + 1 transition. 


(34) 


Wat — l^cav + NqU. 


(35) 
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Figure 3: Selective generation of a TVo = 2 photon (upper panel) and Nq = 3 photon (lower panel) Fock state: Population 7 Tjv as a function 
of N for different pumping parameters. The points are the result of a purely photonic simulation, the lines are a guide to the eye. Left panel 
parameters: for all curves S = —U, 2YLn/Y pU mp = 0.01, and then for each particular curve 2Ti oss /T purnv = 2 10“ 5 (blue solid line), 2 10“ 6 
(green, dashed line), 210 -7 (red, dash-dotted line), 2 10 —8 (magenta, dotted line). 2U/Y pU mp = 10 3 / 2 (blue solid line), 10 2 (green, dashed 
line), 10 5 / 2 (red, dash-dotted line), 10 3 (magenta, dotted line). Right panel parameters: fora 11 curves S = —U, 2Yln/Y pU mp = 0.01, and 
then 2Yi oss /Y purnp = 5 10 —8 (blue solid line), 5 10“ 9 (green, dashed line), 5 10 —10 (red, dash-dotted line), 5 10 —11 (magenta, dotted line). 
2U/Y P ump = 2 10 5 / 2 (blue solid line), 2 10 3 (green, dashed line), 2 10 7 / 2 (red, dash-dotted line), 2 10 4 (magenta, dotted line). The goal of 
these choices of parameters was to control the steady-state ratios P(N + 1)/P(JV) = 10“ 2 and P(N)/P(Q) = 0.1,1,10,100 (blue, green, 
red, magenta). 




Figure 4: Left panel: Purely photonic simulation of the two-time coherence function g ( 2 ) (r) for a strongly nonlinear regime in a (metastable) 
No = 2 photon selection regime. The inset shows a magnified view of the short time region. Parameters: 2 U/Y pU mp = 100, 2Yi oss /Y P ump = 
2 10 -3 , 2Yem/Ypump = 0.2, 5 = — U\ in the language of Fig[3] the present parameters would correspond to a regime where the IV = 0, 2 
states are almost equally occupied. Right panel: Preparation of the metastable state at No = 2 starting from a N = 4 7t( 4) (red dot-dashed) 7r(2) 
(green dashed) 7r(0) (blue solid). Same parameters as in Fig[3] 
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As a result, only this last transition is dominated by emission, while all others are dominated by losses. 

In terms of the diagrams in FigQ] the stationary distribution 7r ; y is therefore sharply peaked at two specific values, 
TV = 0 and at TV = Nq. Examples of this physics are illustrated in Fig[3j the two peaks are always clearly visible, 
but depending on the parameters their relative height can be tuned to different values almost at will. It is however 
important to note that having a sizable stationary population in the TV = TVo peak requires quite extreme values of the 
parameters as population would naturally tend to accumulate at TV = 0 and this difficulty turns out to be exponentially 
harder for larger TVo- 

The physics underlying this behaviour can be easily explained in terms of the asymmetry in the switching mech¬ 
anisms leading from TV = 0 to TV = TVo and viceversa. The former process requires in fact a sequence of several 
unlikely emission events from TV = 0 to TV = TVo — 1 as emission is favoured only in the last step. On the other hand, 
decay from TV = TVo occurs as a consequence a single unlikely loss event from TV = TVo — 1 to TV = TVo — 2: as soon 
as the system is at TV = TVo — 2, it will quickly decay to TV = 0. 

The rate Y acc of such an accident can be estimated as follows: the probability that the system in TV = TVo — 1 
decays to TV = TVo — 2 is a factor (TVo — l)I) 0 ss/(-(Vore m ) smaller than the one of being repumped to TV = TVo- As 
the rate at which the system decays from TV = TVq to TVq — 1 is approximately equal to NqTi oss , one finally obtains 


= N 0 T lo 


(TVo — l)rZoss 
' TVnHL 


<c TV 0 r Zo 


(36) 


This longer time scale r occ = T~^ c is clearly visible in the long tail of the time-dependent that is plotted in the 

left panel of Figj4] The quick feature at very short times corresponds to the emission rate T em . 

If needed, the characteristic time scale r acc could be further enhanced by adding a second atomic species whose 
transition frequency is tuned to quickly and selectively emit photons on the TV — 2 — > TV — 1 transition. In this way, 
the accident rate can be efficiently reduced to T^cc — r; oss (rz oss /rg m ) <C T acc . By repeating the mechanism on k 


transitions, one can suppress the accident rate in a geometrical way to rice — Fz oss (T i oss /r® m ) k -C r acc . Finally, 
the Fock state with Nq photons can be fully stabilized to an infinite lifetime and no problem of metastability if TVo 
different atomic species are included so to cover all transitions from TV = 0 to TV = TVo. 

From a slightly different perspective, we can take advantage of the slow rate of accidents T acc to selectively 
prepare a metastable state with TV = TVo photons even in parameter regimes where the TV = 0 state would be 
statistically favoured at steady-state. Though the state will eventually decay to TV = 0, the lifetime of the metastable 
TV = TVo state can be long enough to be useful for interesting experiments: The idea to prepare the state with Nq 
photons is to inject a larger number TV > TVo of photons into the cavity: the system will quickly decay to the TV = TVo 
state where the system remains trapped with a lifetime T~^ c . 

The efficiency of this idea is illustrated in the right panel of Figj4] where we plot the time evolution of the most 
relevant populations The initially created state with TV = N rn photons quickly decays, so that population 
accumulates into TV = TVo on a time-scale of the order of Fz oss ; the eventual decay of the population towards TV = 0 
will then occur on a much longer time set by Y acc . It is worth noting that this strategy does not require that the 
initial preparation be number-selective: it will work equally well if a wide distribution of N m are generated at the 
beginning, provided a sizable part of the distribution lies at TV > TVo. Furthermore, this idea removes the need 
for extreme parameters such as the ones used in Figj3] to obtain a balance between 7r(TV) and 7r(0): as a result, the 
difficulty of creating a (metastable) state of TVo photons is roughly independent of No- 

These results show the potential of this novel photon number selection scheme to obtain light pulses with novel 
nonclassical properties: for instance, upon a sudden switch-off of the cavity mirrors, one would obtain a wavepacket 
containing an exact number of photons sharing the same wavefunction. With respect to the many other configurations 


discussed in the recent literature to produce TV-photon Fock states and photon bundles (56, 57, 58], our proposal 


has the advantage of giving a deterministic preparation of a TV-photon Fock state in the cavity, which can then be 
manipulated to extract light pulses with the desired quantum properties. 


4. Cavity arrays 

After having unveiled a number of interesting features that occur in the simplest case of a single-cavity, we are now 
in a position to start attacking the far richer many-cavity case. From now on we consider that the isolated photonic 
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Hamiltonian is the Bose-Hubbard one with tunelling J and interaction constant U. Throughout this section, we shall 
make heavy use of the purely photonic description previously derived, which allows to consider bigger systems with a 
higher number of photons. A numerical validation of this approach against the solution of the full atom-cavity master 


equation is presented in Appendix E 


4.1. Markovian regime 

We begin by considering the Markovian limit of the theory, which is recovered for T pump = oo, i.e. for a 
frequency-independent gain. In this case, the emission term of the master equation for photons Eq|TO|reduces to the 
usual Lindblad form 

r o k 

£em = ^2 [ 2a lP a i ~ a i a \P ~ P^a\ . (37) 


i=1 


For a single cavity, the stationary state is immediately obtained as 


ir N = 


( r° \ N 

I em \ 


l-p*- Vr loss) 

k los 


(38) 


a necessary condition for stability for this system is of course that F° m < r/ os , 5 . For r° m > T/ os . s amplification 
would in fact exceed losses and the system display a laser instability: while a correct description of gain saturation 
is beyond the purely photonic theory, the full atom-cavity theory would recover for this model the standard laser 


operation |_52j, [481, 2' 


For larger arrays of k sites, a straightforward calculation shows that in the Markovian limit the stationary matrix 
keeps a structureless form, 

Poo = ^ TTjV^iV, (39) 


N 


with 


7T/V = 


Em d m (pEj) 


M 


N 


(40) 


Here, Dm = ' s t ^ le di mens i°n of the Hilbert subspace with a total number of photons equal to N and ly is 

the projector over this subspace. The interested reader can find the details of the derivation in |Appendix C| 

This result shows that independently of the number of cavities and the details of the Hamiltonian, in the Marko¬ 
vian limit the density matrix in the stationary state corresponds to an effective Grand-Canonical ensemble at infinite 
temperature (3 = 0 with a fugacity z = = T^ m /Yi oss determined by the pumping and loss conditions only: All 

states are equally populated and the system does not display much interesting physics. In particular, the steady state 
does not depend on the tunelling amplitude J and on the photon-photon interaction constant U 


4.2. Effective Grand-Canonical distribution in a weakly non-Markovian and secular regime 

The situation changes as soon as some non-Markovianity is included in the model. In this section we start from 
a weakly non-Markovian case where all relevant transitions adding one photon have a narrow distribution around the 
bare cavity frequency, | uff — uj cav | ■C F pum p- We also assume a secular limit where U, J 3> r° m , Fj oss , so that 
the non-diagonal terms of the density matrix in the photonic hamiltonian eigenbasis oscillate at a fast rate and are thus 
effectively decoupled from the (slowly varying) populations. In this limit, we can safely assume that all coherences 
vanish and we can restrict our attention to the populations. This somehow critical approximation will be justified 
a posteriori in the next section, where we treat perturbatively the coupling of populations to coherences and show 
both analytically and numerically that in the weakly markovian regime, their contribution is of higher order in the 
’non-markovianity’ parameter 1/T pump and therefore can be safely neglected. 

Under these assumptions, the transferrate on the |/') — > \ f) transition where one photon is lost from N + 1 to N 
has a frequency-independent form 

Tf-yf = Floss \(f\ a\f)f , (41) 
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while the reverse emission process depends on the detunings A ff = ojff — u) cav and 5 = oj cav ~ oj a t as 


?>-►/' = r" ro I </Vl/> 


Wcav W a t)^ 


ri, 


^r“ m |(/| a t|/')r i-pA f , f + o(A f , f y 


( 42 ) 


with 


(Tpump/zy 


{.Wcav Wat “1“ (r pump /'2) 

2(uj 

cav Wat) 


( W C av tU at “1” {J'pump/ty 


2 • 


(43) 

(44) 


In this expression, the weakly non-Markovianregime is characterized by having \3Af/f\ <C 1: in this case, the square 
bracket in Eq|42lcan be replaced with no loss of accuracy by an exponential 

1 — /3A ff ^ e~P A ft, (45) 


which immediately leads to a Grand-Canonical form of the stationary density matrix 


Poo = 


(46) 


with an effective chemical potential 


1 / f 0 


■ W c 


(47) 


and an effective temperature fc^T = 1 /3: most remarkably, even if each transition involves a small deviation from 
the bare cavity frequency u> cav , the cumulative effect of many such deviations can have important consequences for 
large photon numbers, so to make the stationary distribution strongly non-trivial. Remarkably, both positive and 
negative temperature configurations can be obtained from Ea l44li ust by tuning the peak emission frequency uj at either 
below or above the bare cavity frequency uj cav . As expected for a thermal-like distribution, detailed balance between 
eigenstates is satisfied 


= \{f'\ ^ |/) 


-0w f , + 



= 0, (48) 


but it is crucial to keep in mind that this thermal-like distribution does not arises from any real thermalization process, 
but is a consequence of the specific form chosen for the pumping and dissipation. The application of this concept to the 
study of effective thermalization effects in a driven-dissipative non-Markovian condensate in the weakly interacting 
regime will be the subject of a future work, also with an eye to photon [34] and polariton fi~4l [33 1 Bose-Einstein 
condensation experiments. 

A numerical test of this result for a two cavity system with a strong pumping r pump U■ ■) and a large enough 
photon number so to induce appreciable nonlinear effects is shown in FigJ3 The results of this comparison are 
displayed in the left and central panels: excellent agreement between an exact resolution of the photonic master 
equation and the grand canonical ensemble ansatz is found in both the average photon number and the first-order 
coherence. 
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Figure 5: Left and center panels: average number of photons n\ = {a\a\) (left) and spatial coherence g^\ = { a \ a 2 ) / {o\a-\) (center) in a 
two cavity system with small U /T purnp and J/T purnp as a function of the non linearity U at fixed T purnp ■ In red dots, exact resolution of the 
photonic master equation, and in black solid line the grand canonical ensemble ansatz. Parameters : 2J/T pU mp = 0.02, 2r / oss /Tp Ump = 0.002, 
2T era/Tpump = 0.0014, 25/T purnp = 0.6. Right panel: purely photonic simulation of the relative quantum coherence between two arbitrarily 
chosen two-photon eigenstates pij / y/ pupjj as a function of 1 /T pU m P (the result does not depend on the specific eigenstates considered). As 
expected, this coherence vanishes in l/Fp Ump in the Markovian limit 1/r pU mp —> 0. The value above 1 for large 1/T purn p signals breakdown 
of positivity of the density matrix as we move out of the validity regime of the purely photonic master equation. Parameters: J/r Zoss = 1, 

r e m/r Zoss = o.5, s = — r Zoss , f//r Zoss = 2 . 


4.3. Beyond the secular approximation 

In the weakly non-Markovian regime, the validity of the effective Grand-Canonical description can be extended 
outside the secular approximation according to the following arguments. As a first step, we decompose the master 
equation as 

^ = [M 0 + SM]p, (49) 

where the super-operators Ml and 6Ml act of the linear space of density matrices p as 


M 0 [p\ = -i [ H , p\ + [ 2a i 


pa f - a\atp - pa\at 


2=1 


r u r 

1 em I A f + ^ A f 

— 2 ^ [ a iP a * + a iP ai ~ a i a iP~ P a i a l 


i=1 


and 


with 


and 


»=1 


r ° K r 

5Ml[p] = —+ a\p5ai - atSajp - pSa l a\ 

(«! + <V) i 


F° 

~t _ em 


r° 


'i ^ i a - _VM /f\j 




(/VII/) 


from which we deduce that 


/ pump—too \ 1 pump \ -L pump J 


, (50) 


(51) 


(52) 


(53) 


(54) 


Using similar arguments to the Markovian case of | Appendix C| we can easily show that the grand canonical distribu¬ 
tion is a steady state of this modified Ml o operator. 


M 0 {e 


@Nti —pH \ _ 


) = 0 , 


(55) 
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As the correction term 8A4 vanishes in the Markovian limit proportionally to l/T pump , we can calculate the lowest 
order correction to the steady state in 8A4. Expanding the steady state in powers of 1/T pump keeping a constant 
(w co „ — oj a t )/r p U mp, we see easily that the first order corrections in eq. (l54l > are purely imaginary so that populations 
are perturbed only to second order in BA p f. In our Markovian limit, these corrections then vanish even if we perform 
simultaneously the Markovian and thermodynamic limit. 

Secondly, coherences (which are exactly zero in the Markovian case, see Sec 14. II ) should be then proportional to 
1/r p U mp- However, we have shown in Appendix D that the linear contribution to coherences vanishes when we sum 
over all sites of the system. We conclude thus that in the weakly non Markovian limit, coherences between eigenstates 
of the hamiltonian are quadratic in 1 /T pump and therefore remain very small even out of the secular approximation. 

As a further verification of this analytical argument, in the right panel of FigJJwe have shown the T pump de¬ 
pendence of the coherence between an arbitrary pair of two-photon states as well as the error in the population of 
an arbitrary eigenstate, between the true steady state and the grand canonical distribution. As expected on analytical 
grounds, both these quantities scale indeed as 

From these arguments, we conclude that the breakdown of the secular approximation which occurs in the thermo¬ 
dynamic limit where the spectrum become continuous should not affect the effective thermalization of the steady state 
in the weakly non-Markovian regime of large T purnp . Even if the steady-state is not affected, we however expect that 
the relatively strong dissipation will significantly affect the the system dynamics. A complete study of this physics 
will be the subject of a future work. 


5. Two cavities with strong non linearity 


5.1. Towards Mott-insulator physics 

As a final example of application of our concepts, in this last section we present some preliminary results on the 
most interesting case of two strongly nonlinear cavities with U A> Tpump- extending the photon-number selectivity 
idea to the many-cavity case, we look for many-body states that resemble a Mott insulator 117 . H |49j. As in the 
single cavity case, the strong pumping T em 3> Ti oss would favour a large occupations of sites, but is counteracted by 
the effect of the nonlinearity U I pump which sets an upper bound to the occupation: the result is a steady-state 
with a well-defined number of photons per cavity. 

The result of numerical calculations based on the photonic master equation are shown as black lines in FigJ^a-c) 
in the cc r:av = ui a t case: for a high emission rate T° em and a strong non linearity U, signatures of the desired Mott 
state with one particle per site are visible in the steady-state average number of photons that tends to 1 for a strong 
nonlinearity U [panel (a)], in the probability of double occupancy that tends to 0 [panel (b)], and in the one-body 
coherence between the two sites that also tends to 0 [panel (c)]. 

While these results are a strong evidence of Nq = 1 Mott state, a similar calculation for larger No > 2 Mott states 
is made much more difficult by metastability issues and the Mott state would typically have a finite lifetime. As in the 
single cavity case, we expect that this problem could be fixed by adding several atomic species on resonance with the 
different photonic transitions below Nq. 

Based on this preliminary analysis, we can attempt to make some claims on the structure of the non-equilibrium 
phase diagram of our model. As for J = 0 one can efficiently create a Fock state in each cavity, we expect that for 
small J the system will remain in a sort of Mott state. On the other hand, in the weakly interacting regime we expect 
the system to display a coherent Bose-Einstein condensate |50i|. In between, one can anticipate that system should 
display some form of non-equilibrium Mott-Superfluid transition. Analytical and numerical studies in this direction 
are in progress. 


5.2. An unexpected mechanism for coherence 

The red dashed lines in the same panels Fig[6[a-c) show the same simulation for a weaker emission rate r° m , which 
allows to consider weaker values of the nonlinearity without increasing too much the photon number. In particular, in 
panel (c) we see that the non-negligible value of 2J/T pump is responsible for a significant spatial coherence between 
the two sites, which attains a maximum value g[ J, ss 0.26 for an interaction strength 2(7 /T pump ~ 0.16 of the same 
order of magnitude as the tunnel coupling 2J/T pump = 0.2. 
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Figure 6: Purely photonic simulations of steady-state observables as a function of 2U/T pU mp in a two-cavity system: (a) average number of 
photons ni = ( a[a ±), (b) one-site two-body correlation function g^\ = (ajajaiai) = (ni(ni — 1)), (c) inter-site one-body correlation 
function = (a\a 2 )/(a[ai). Parameters: 2 J/Y pU mp = 0.2, 2Ti oss /Y pU mp = 0.002, 2T e m/^'pump = 0.06 (solid black line). Red 
dashed line, same simulation with a weaker 2r ern /T purnp = 0.00144. Panel (d), from left to right: state occupancy, energy and two site spatial 
coherence of the different eigenstates of the hamiltonian, at the maximum coherence point 2 U/T pU m P = 0.16 of the red dashed line. 
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The quite unexpected appearance of this coherence can be understood as follows. On one hand, in the absence 
of tunneling J = 0, all the dynamics is local and we do not expect any spatial coherence. On the other hand, in 
the absence of interactions U = 0 and for zero detuning, symmetric and anti-symmetric states are equally close to 
resonance (albeit with opposite detuning) and then equally populated, so there should not be any coherence either. 
However in presence of both tunnelling and small interactions (i.e. for J,U ^ 0 and U <C J), the energy of all 
eigenstates (symmetric/anti-symmetric states with various photon numbers) is perturbatively shifted in the upward 
direction by (small) interactions U. As a result, symmetric states, which are below the resonance, get closer to 
resonance and become more populated than the anti-symmetric ones, which get farther to the resonance and are thus 
depleted. As one can see in the plot of the energy, the spatial coherence and the steady-state occupancy of the different 
eigenstates shown in Figj6jd) for the maximum coherence point, this induces an overall positive coherence between 
the two sites. 

Even though the nonlinearity is only active for states with at least two photons, it is interesting to note that also 
in the N = 1 manifold the antisymmetric state is less populated than the symmetric one. This population unbalance 
is inherited from the one in the above-lying N > 1 states, as the decay preferentially occurs into the symmetric state. 
Since no coherence is expected in both limiting cases of purely interacting U J and non interacting U /0 = 0 
photons, the maximum of the coherence is obtained when interactions and tunelling are of the same magnitude, 
U w J: this result is clearly visible in panel Figjhjc). 

Investigation of this many-body physics in the more interesting case of larger arrays which can accommodate a 
larger number of photons requires sophisticated numerical techniques to deal with the dynamics in a huge Hilbert 
space [£9j, 60] and will be the subject of future work. A very exciting advance in this direction was recently published 
in f38tl for strongly interacting photons in the presence of a synthetic gauge field for light: analogously to the Mott 
insulator state studied here, the combination of the effectively frequency-dependent pumping (obtained via a two- 
photon pumping in the presence of an auxiliary lattice) and the many-body energy gap was predicted to generate and 
stabilize fractional quantum Hall states of light. 


6. Conclusions 

In this work we have proposed and characterized a novel scheme to generate strongly correlated states of light in 
strongly nonlinear cavity arrays. Photons are incoherently injected in the cavities using population-inverted two-level 
atoms, which preferentially emit photons around their resonance frequency. The resulting frequency-dependence 
of the gain will be the key element to generate and stabilize the desired quantum state. A manageable theoretical 
description of the system is obtained using projective methods, which allow to eliminate the atomic degrees of freedom 
and describe the non-Markovian photonic dynamics in terms a generalized master equation. 

The efficiency of the our pumping scheme to generate specific quantum states is first validated on a single-cavity 
system: for weak nonlinearities, a novel mechanism for optical bistability is found. For strong nonlinearities, Fock 
states with a well-defined photon number can be generated with small number fluctuations. 

In the general many cavity case, in the weakly non-Markovian case the steady-state of the system recovers a 
Grand-Canonical distribution with an effective chemical potential determined by the pumping strength and an effective 
inverse temperature proportional to the non-Markovianity: This very general results may have application to explain 
apparent thermalization in recent photon and polariton condensation experiments. 

Finally, the power of a frequency-dependent pumping to generate strongly correlated states of light is illustrated 
in the case of a strongly nonlinear two-cavity system which, in the strongly non-Markovian regime, can be driven 
into a state that closely reminds a Mott-insulator state. A general study of the potential and of the limitations of 
the frequency-dependent gain to generate generic strongly correlated states with many photons will be the subject of 
future work. 
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Appendix A. Derivation of the purely photonic master equation via projective methods 


In this Appendix, we give more details on the derivation of the photonic master equation ®. Starting from the full 
atom-cavity master equation ©, we show how for a sufficiently small atom-cavity coupling SI r the atomic degrees of 
freedom can be eliminated. The frequency-dependence of the atomic amplification is then accounted for as a modified 
Lindblad term (ITOl) . Our treatment is based on the discussion in the textbook 14811. 


Appendix A. 1. General formalism 

We consider a quantum system which undergoes dissipative processes. As it is not isolated, its state can not be 
described by a wave function but by a density matrix p evolving according to the master equation : 

d t p = C(p(t)), (A.l) 

where C is some linear “super-operator” acting on the space of density matrices. Given an arbitrary initial density 
matrix p(to), the density matrix p at generic time t is equal to p(t) = e c ^~ to ’ p(to) ■ 


Now we are only interested in some part of the density matrix, which can represent some subsystem. This can be 
described by a projection operation on the density matrix Vp . We call Q = 1 V the complementary projector. We 
decompose the Lindblad operator C in two parts Co and 8C such that: 

r C = C 0 + 5C 

< VC q Q=Q£oV = 0 (A.2) 

V6£V = 0. 


Such a decomposition is always possible. 

Then we define a generalised interaction picture for the density matrix and for generic superoperators A with 
respect to the evolution described by the free Co and the initial time to'- 


p(t) = e to ^p(t) 

A(t) = e - £ o(*-to)_4 e £o(t-to). 


(A.3) 


As discussed in 14811 , we can get an exact closed master equation for the projected density matrix in the interaction 
picture 


which translates to 


dtPp{t) = [ dt'£{t,t')Vp{t'), 

Jt 0 

dtVp(t) = Co(p(t))+[ dt'£(t-t')Vp(t') 

J to 


in the Schrodinger picture. In the interaction picture, the self energy operator E is defined as: 


00 nt rt\ pt n — l 

s (m') = E/ / ■■/ 

•v, _o " t f J t' J t' 


dti..dt n V6C(t)Q8C(ti)Q6C(t 2 )...Q6C(t n )Q5C(t')V 


(A.4) 


(A.5) 


(A.6) 


and results from the coherent sum over the processes leaving from V, remaining in Q and then coming back finally to 
V. In the Schrodinger representation, we have : 


E(t - t') = e £o( *“‘ o) E(f, = S (o, t' - t)e L °( t ~ t ''>. 


(A.7) 


We call t c = 1/A u) the characteristic decay time / inverse linewidth for the self energy, which corresponds 
in general to the correlation time of the bath, and we estimate the rate of dissipative processes as T ~ Er c = 
f t dtE(t. to)- We put ourselves in the regimes in which, with respect to these dissipative processes, the bath has a 
short memory, ie T <C A lj. In that regime the density matrix in the interaction picture is almost constant over that 


18 




time t c . Furthemore, if t — to t c then the integral in eq (lA.4b can be extended from — oo to t. From this equation 
and from (1A.4I >. we get an equation of evolution for the density matrix which is local in time : 


dtPp[t) 


pOO 

/ c?rD(t, t — r)Vp(t) 

Jo 

e -£(i-to) S ( 0 ) _ T ) e £(i-to) 
poo 

/ drD(0, — r)Vp(t). 

Jo 


L 


dr 


o -C(t-t 0 ) 


In the Schrodinger picture this gives the time-local master equation : 


(A.8) 


d t Vp(t) 


Co 


f 


drS(0, — t) 


Vp{t) = C eff Vp{t), 


(A.9) 


with 

pOO 

C ef f =£ 0 + drl](0, -r). (A. 10) 

Jo 

It is worth stressing that while the bath is Markovian with respect to dissipative processes induced by the perturbation 
/q 30 drE(0,—r) , no Markovian approximation has been made with respect to the dynamics due to jOq, which can 
still be fast. For the specific system under consideration in this work, this means that the emission rate r em has to 
be slow with respect to the gain bandwidth set by the atomic pumping rate T pump , which is the case in the weak 
coupling limit y/N at fln -C F pump , but no restriction is to be imposed on the parameters U, J and uj cav — 0J a t of the 
Hamiltonian, which can be arbitrarily large. This means that the physics can be strongly non-markovian with respect 
to the Hamiltonian photonic dynamics. 


Appendix A. 2. Application to the array of cavities 
Preliminary calculations 

With the notation from section^ we choose the projectors in the form : 


Vp 


3 (i)„(i)„(i) 

'i 


,(i)„(i)„(i) 

'i 


1 Tr at (p ), 


(A.11) 


where we have performed a partial trace over the atoms, and then make the tensor product of the density matrix and 
the atomic density matrix with all atoms in the excited state. We chose this particular projector because in the weak 
atom-cavity coupling regime, we expect atoms to be repumped almost immediately after having emitted a photon in 
the cavity array, and thus to be most of the time in the excited state. Moreover this projection operation gives us 
direct access to the photonic density matrix, and thus we do not lose any information on photonic statistics. With the 
notation of the previous section we have : 



£(p) — —i [H p h + Hat + Hi, p] + Cdi S s(p\ 

(A. 12) 

with 

Cdiss Cp Urn p, a t Class,cav 

(A.13) 

We decompose C in two contributions. The first one is : 



£-0 (p) — [Hph + Hat J p\ + £loss,cav{p ) — -4(p) + V AQ(p) 

(A. 14) 

with 

-p k Nat 

, LpumpN-^N-' \ -( l ) +(l) -( l ) +(Z) 

AP)= 2 [ a i a i P + P a i a i ■ 

(A.15) 


i=l 1=1 


19 










The superoperator Co verifies the condition (lA.-l i: The last term in the expression of eci. dA. 1 41 ) comes from the fact 
that the pumping term A in Co does not verify this condition: as a result, we have to remove the part unfixed by 
projector and put it in the other operator : 


p k Nat 

SC(p) = i\H I ,p\J r 2<j+ (l) pa i {l) - VAQ(p). (A.16) 

i=l Z=1 

These two operators then satisfy to the conditons (1A.2I) . and we can apply the projection method to get the evolution 
of Vp(t), that is of Tr at (p)(t). As we are interested in the regime in which T pump y/NatflR, Tj oss , we will 
compute the self energy at the lowest non zero order of these two latter parameters. Since l’/ os ., quantifies the photonic 
loss rate, we will approximate the photonic dynamics as being a Hamiltonian one during the time while the atom is 
reinjected in the excited state, ie during the characteristic time \/T pump of the integration kernel of eci. dA. 5I >. To this 
order of precision, the calculation for one cavity is easily generalizable to k cavities, thus we will restrict for simplicity 
to the case of a single cavity containing a single two-level atom, N at = 1. 


Self energy calculation : 

We are going to calculate the self energy to the lowest order in f Ir. We have 

SC = £ pump - i{H+ + H~) l + i(H+ + H~) r - VAQ , (A.17) 


with 


*~~ / pump 

H+ = 
H~ = 


(p) — Cpump&+ P& 
ClR<j + a 
flRa~a t 


(A.18) 


By (H^r/r we intend the superoperator multiplying a matrix p by the matrix H ± on its left/right. First we have 
Cpump'P = VAQV = HrV = HffP = 0, so starting from a projected state V p, we have to start with Hf or H^. 
In fact to the lowest order in tf,R the non zero contributions to the self energy are : 


A = -t)V 

B = -VH^H+lt' -t)V 
C = VH+Hf{t' ~t)V 
D = VHfH+it' - t)V 

E = $1 dtVC purnp (t)QH+{t - - t)V (A.19) 

F = Si dtVC pump (t)QH+(t - t)Hn(t' - t)V 
G = - Si dtVAQH+(t- t)Hf(t' - t)V 
H = — Si dtVAQHf{i-t)H+(t' - t)V, 

with 

£(0, t' -t) = A + B + C + D + E + F + G + H. (A.20) 

We then calculate the different processes, applied on some projected matrix Vp: 
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A(Vp) = 
B{V P ) = 
C(Vp) = 
D(V P ) = 

E(Vp) = 


-n 2 R e^ at ~ Tpump,2){ - t ~ t ' ) aa ] (t' - t)Vp 
--Uf ^ e _(iWa ‘ +^p “ mp/2)(^_^,) 7 ? pa(i , - t)a t 
Ct R e^- iUat _r !> u "*p/ 2 ) (*-*') at (t’ - t)Vpa 
fl R e^ iUat+I ' pump ^ 2 ^ t ~ t '^cJVpa(t' - t) 


r o 2 

x pump* L R 


f 


dt e( r pump / 2 )(t t ) 


e(“ a * r P“ m p/ 2 )(* a){t'- t)Vpa(t - t) 


f{v p ) = r 


O' 

pump * L R 


■L 


dt e 


(iuiat-F P ump/2)(t—i) 


g( iUai r pump / 2 )(t t') a + [t - t)Vp a {t'- t) 


G(Vp) = 


H(Vp) = 


—V o 2 

x pump * L R 


dt i ,UJat ^pump /2) (t t) 


r pump / 2 )(t £ )/ 7 .lY/ / _ 

rt 


eK .~az zpump,-,^ - -t)Vpa{t -t) = -E(Vp) 


__ p n2 
X pump * L R 


L 


die {iUai 


-r„ 


p/ 2 )(t-t) 


e ( r p „ mp / 2 )(i t') a t(£ _ t)Vpa,(t' - t) = -F{Vp) 


(A.21) 


where by a[t' — t ) we intend the evolution of the photonic annihilation operator in the photonic hamiltonian interaction 
picture (we remind that we neglected photonic losses during the integration time). We see that the last four contribution 
cancel each other, and that only the first four contributions remain. 


Master equation 

Using the expression for the self-energy £(f) derived in the last section, as well as general results on the master 
equation obtained by projective methods in Sec jAppendix A.l we then obtain the (temporally non-local) master 
equation : 


d t Vp = -i [H ph ,Vp\ +C T {Vp) 

+ n% f d r e <—™/»V - r>) a(-r) 

+ft 2 R j dr e -(i«.*+r,W2)T 0 t(_ r ) (e Ca ^Vp{t - t)) a 

— f l 2 R J dr e^‘ Uat ~ rpv - rnp ^ 2 ^ T aa^(— t) (e c °^Vp(t — r)^ 

-fl R f dr e~ < - iuJat+I ' pump/2 ' >T (e^^Vpit - t)^ a(-r)a^. 


(A.22) 


21 





At lowest order in Qr, we can assume the interaction picture density matrix in the convolution product to be constant, 
p(t — t) ~ p(t), i.e. e c ° T p(t — t) ~ p(t). Making the trace over the bath we get: 

df.Pph — ^ [Hphi Pph\ A t-'T^Pph) (A.23) 

POO 

+fl 2 R / dr e ( ' iUat ~ Vpnmp,2 ' )T a) (~T)p ph {t)a 

J 0 

pOO 

+n 2 R / dr e~ {iuiot+r ^ /2)T a^P ph {t)a{-T) 

J 0 

POO 

-fi 2 j / dr e { - iuiat ~ Tpump/2)T aa) (-r)p ph (t) 

Jo 

pO O 

/ dre“ (i ““ t+rp “ mj)/2)r p p / l (<)a(-r)a t , 

io 

then we can perform completely the integral and we get our final form for the non Markovian master equation, which 
is local in time : 


OfP — ^ \Hphi Pph\ 


[2apa^ — 


a) ap 


pa) a] 


+ 


2 n 2 R 


- pump 


[a) pa + a) pa — aa) p — paa)\ 


(A.24) 


with 



J 0 °° rf T e(-”“'“ r p«™p/ 2 ) T a(_T), 
rppmp / 0 °° dr e^ lUat ~ Tpump / 2 ' )T a < {— r) 



(A.25) 


where a(—r) means the photonic annihilation operator in the photonic hamiltonian interaction picture. 

If |/) and |// are two eigenstates of the photonic hamiltonian with a photon number difference of one, we see 
that the matrix elements of the modified annihilation and creation operators a and a) involved in the emission process 
are : 


(/I ^ m 
</'|a|/> = 


_ r pump / 2 _ 

i(uJat ^ f f f pump 

_r 'pump /“^ _ 

^ ff' )+r pump / 2 


72 (/I at in 
(/I a |/). 


(A.26) 


The non-Markovianity comes from the energy-dependence of the prefactors. 

For several cavities the reasoning is exactly the same and we get the multicavity master equation : 


dtp = -i [H ph , Pph ] + yp [ : 


2atpa\ - a\a t p - pa\at 


i-1 


2 n 2 R 


n 

i=l 


a\pai + a] pat - ctidjp - pdta\ 


(A.27) 


with 


(/I a, m 

(/'|aj |/) 


pump 


/ 2 


i'iyiat ^/'/) “I” r pump /2 


(/I Oil/') 


pump 


/2 


'i(y)at ^/'/) “1“ r pump /2 


(/'I at I/). 


(A.28) 
(A.29) 


where here also |/) and // are two eigenstates of the many cavity photonic hamiltonian: once again the emission 
depends on the many body photonic dynamics via the prefactors in (lA.28HA.29l l. 


Appendix B. Lindblad form for the photonic master equation in the secular approximation 

In this Appendix, we present the derivation of the Lindblad form Eq.[l7]for the photonic master equation including 
non markovian effects under the secular approximation. To do this, we calculate the matrix elements C em j, p ^ j of 

the emission superoperator coupling the term of the density matrix in the eigenstate basis (/1 p\f) to (/' | p\f), under 
the assumption Aw = w^, p — w^ j ~ 0, as explained in Sec. 12.31 
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Calculation of the a\pa,i + a\pat contribution: 


(/'I 5t |/) (/| Pl/X/M/') + (/'I ^ |/) (/| p|/)(/|a i |/ / ) 

r/| I A / f \ „| ;\/>UJ h\ I Fpump/Z 


= <mi/> </M/X/H/'> 


^ pump/^ 


“I” pump/% i'iyiat ^ pump /2 


• (B.l) 


Considering that under the approximation Aw ~ 0, we have that - i{uiat -^)+ Vpump /2 ~ 
we obtain thus the following contribution: 

(/'I a! I/) (/I pI/X/N/') + (/'II/) (/I p|/X/N/'> 

r/l t | f \ _ r p«mp /2 


— (/I a I I/) 


(/|P|/)- 


rpump/2 


X/M/'> 


W /',/) 2 + (r pU mp/ 2) 2 y/ (Uat ~ u f'j ) 2 + (r p ump/ 2) 2 

= (/ , |4|/)(/|p|/)(/|a i |/'), (B.2) 

with a,; defined in Eq{l9] We see that the "imaginary" contribution cancels out, and that the "real" contribution has 
been divided in two multiplicative contributions on the left and the right of the density matrix. 

Calculation of the ajajp + pa^aj contribution 


Let us calculate the left product: 

(/'I oil/") (/"I aj |/) </|p|/) = (/'k|/">- 
= </VX/"> 


. -To (/"I “I I/) </IpI/> 

i ^ pump/^ 


(r pump /2 ) 2 


— Z - 


M at)^/pump / 2 


</"l4 I/) </IpI/)- (B.3) 


_iftat ^/",/X d” (B p um P /2) 2 ^/",/X -f (r p1imp /2) 2 _ 

Considering that under the approximation w/'j ~ 0, we have that uJf" j — and so: 


(Tpump/2) 2 


(Tpump/2) 2 


(B.4) 


[ptat “f (Bpump/2) y/( UJ a t ^f",f') 2 “h (r p um P /2) 2 \f{ptat ^f",f') 2 “h (B p itm P /2) 2 

As a consequence: 

(/'I Oi I/") (/"I «ll/) </IpI/> 

* -i (/'I Oi ID (^/C/-^Xr„/2 (r| a t |/} (/ | pl ~ f) + {f ,| s . !n {r | 5 t | /} {/ | p |/). (B.5) 

V^aZ '\*-pump/^) 

Finally, let us calculate the right product: 

r pump/^ 


{f\p\f){f\a>i\f"){f"\al\f') = (f\p\f)~ 
= (f\p\f) 


B pump /2 


</l«i|/"X/ ,, l«II/ / > 


(r pwmp /2) 2 


+ Z \^/ // ?/ ^at)rpnmp/2 

(k^aZ (k^aZ ^/ /; ,/)^ (Bpnmp/2) 2 

As before, coj„ j cop, p, so: 


</|ax/"X/ ,/ l«II/ / > (B.6) 


(r pump /2) 2 


(r pump /2) 2 


(w ot w/» ;) 2 + (r pump /2) 2 ^(wat - w^,^) 2 + (r pump /2) 2 ^/(w ot - wf,,j,) 2 + (r pump /2) 2 ’ 

j"f tu a ^)r pump /2 {top, w a t)r pump /2 

Wj// ^) 2 + (r pump /2) 2 Wj - // y/) 2 + (r pump /2) 2 
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(B.7) 































and thus 


(f\p\f)(mh(f"\al\r) * (/i p\f)(f\°‘i\f"){f , '\®t\f > ) 


+i(f\p\f)(f\ai\n 



( f"\4\f’)■ (B.8) 


(tUat "i~ (Bpump/2)^ 


pump 


Here again, the real part has been divided in two multiplicative contributions, and the imaginary part has been swiped 
to the creation operator on the right. So wether we consider the contribution acting on the left or on the right of 
the density matrix, the imaginary contribution is always carried by the creation operator, so the density matrix is 
multiplied by the same operator on the right and the left up to a minus sign, which gives an anticommutator and thus 
an hamiltonian contribution due to the Lamb shift. 

Sum of the various contributions 

To summarize, keeping only relevant transitions we can consider that the emission dynamics is equivalent to a 
contribution —i [JT Pph\ + 4m in the master equation, with 



(B.9) 



(B.10) 


)</>!!/>, 


(B.ll) 


/' 


which demonstrates the statements of Sec l2.3l 

Appendix C. Exact stationary solution for Markovian case 

In this Appendix, we present a proof of our statements in Sec 14. II We are looking for the steady state for the 
Markovian quantum dynamical process : 


d t p = -i [H, p(t)} + Ci oss + C, 


(C.l) 


with standard Lindblad operators : 



(C.2) 



(C.3) 


i=l 

We want to demonstrate that the following density matrix is an exact steady state : 



(C.4) 


N 


with 



(C.5) 
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First, since the hamiltonian preserves the total photon number, and that the density matrix is equal to the identity 
on each sub-space with a defined photon number, we get that [H, p^} = 0. Second, for the Lindblad operators the 
non-hermitian hamiltonian terms have a simple action on the density matrix : 


pooE a E = 

i 

Poo^ — Npoo — ^ ^ O j iPoo’i 
i 

(C.6) 

Poo ^ ^ O'iO'i 

Poo(N + /c) — ^ ' CLi Cl i Poo ? 

(C.7) 

i 

=aj ai + 1 

= (N+fc)p=o 



where k is the number of cavities. We are left with the special terms of the form pa and apa\ for which we find 
that: 


and 



EE E i/)(/i-(/i4i/') Eipooi n (/ioil/) 

i N n eq ( N -i) 6 f'J' 

ff(N-l) 

EE E I/) (/I' 7r JV-i</|Oi \ f) (f \ a-i |/) 

* N 

f'(N-l) 

E E i/x/i-^-iaiE^i/) 

" /, f(N) > -u-, 

=NfS fJ , 

E E i/> </i ■ 

N f(N) 


E a iP°° a l = EE( ]V + 1 + k)n N+1 |/) (/| . 

i N f(N) 


If we sum all contributions together, it is immediate to see that we get a total zero contribution : 


(C.8) 


(C.9) 


k^loss^Poo) T ^-em(Poo) — 


= EE i/x/i (E e 

N f(N) \ 


^n-i ~ NTiossKn + (N + k)Ti oss 7r N+ i - (N + k)T em n N = 0, (C.10) 


=o 


=0 


which proves our statement. 


Appendix D. Perturbative corrections to the coherences in the weakly non Markovian regime 


In this Appendix we show that the lowest-order correction to the coherences between eigenstates (null in the Grand 
Canonical ensemble of Sec l4.2b are quadratic in the inverse pumping rate T~^ mp and not linear as a naive pertubative 
expansion would suggest. To this purpose, we calculate the first order contributions to the coherences of the operator 
SA4 [defined in eas. lBTl) and (l54t | applied to the grand canonic density matrix and show them to be 0. Let us calculate 
first the contribution of the first two terms : 


E(/i^W i |/ , >= E (f\ Sa l\f)(f'\<h i/') x (/ipooi/o 
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Figure E.7: Comparison of the analytical prediction of the photonic theory (solid black line) to the numerical solution of the full atom-cavity master 
equation (open red points). Stationary value of the average number of photons as a function of the photon loss rate rj oss (left) and of the atom- 
cavity coupling (right). Parameters . 2U /T'pump — 2, 2^lj^/T'pump = 0.02 [left panel (a)], 2U /T pump — 0.6, 2T loss /T pump —— 0.02 
[right panel (b)]. In all panels, 28/T pU mp = 8. 


In the same way : 


Then we know that 


Y (/I a \PooSai \f) = Y (/I a l\f)Cf\^i I/') (/|poo|/) 
i ij 


(f\Se\\f) = — (f\at\f)+0 


(D.2) 


(D.3) 


Let us choose a reference state |/ 0 ) with the same photon number as |/). Then (f\poo\f) = {fo\ Poo \fo) +0(T p ^ mp ). 
All these additional terms give second order contributions, and we do not consider them. Thus to the first order : 


Ei (f\ Sa iP< 


a* + ajpooScii I/'} 

= Y.iJ {f\ a \\ f)0\ a i I/') (MPoo |/o) 


= A 
= A 
= A 


r pump 
r pump 
r pump 


Ei,/(/l4l/)(/l a i I/') 

E* (f\a\ai\n 

(.f\N\f) 


=NfS t f> 


= A 
= 0. 


f -N f S 


f°ff 


(D.4) 


A similar reasoning allows to show that 

Y (/I aiSajpoo + PoeSa-id] \f) = 0 (D.5) 

% 

which completes our proof. 


Appendix E. Further numerical validation of the photonic master equation 

One cavity case 

Here we compare the analytical prediction for the stationary state of the atom-cavity system discussed in Sec [3] to 
a numerical solution of the full master equation Eq.(|4]). For example, in the left panel of Fig lE.7l the stationary value 
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Figure E.8: Comparison of the analytical prediction of the photonic theory (solid black line) to the numerical solution of the full atom-cavity 
master equation (open red points) for a two-cavity system. Stationary value of the average number of photons in the first cavity as a function of 
the photon loss rate F/,,,, (left) and of the atom-cavity coupling f 1r (right). Parameters: 2U/T pllrn p = 7, 2(1 R/T pump = 0.02, (left a) panel); 
2U/T'pump —— 28, /T'pump —— 0.002 (right l>) panel). In all panels, 2.//F pump — 4 and 5 — 0. 

for the average photon number is plotted as a function of the photon loss rate 1 - As expected, the purely photonic 

approach based on the projective method gives very accurate results as long as the pump rate T pump (i.e. the inverse 
autocorrelation time of the atomic bath) is much faster than the loss rate l'/ 0 . s . s . 

A similar plot of the average photon number as a function of the atom-cavity coupling f 1r is shown in the right 
panel. Outside the small Qr regime, the photonic theory tends to overestimate the photon number. This deviation can 
be explained as the theory assumes the atoms to be always in their excited state ready for emission and neglects the 
possibility of an atom reabsorbing the emitted photon before being repumped to the excited state. 

Two cavity case 

Here we give further validation to the purely photonic description used in Sec0]by comparing its predictions with 
the numerical results for the full atom-cavity master equation in a two cavity case. An example is shown in Fig lE.8l 
as in the single cavity case, the agreement is excellent at large T pump and gets deteriorated when T pump is decreased 
to values comparable to L’/„, s . s [panel (a)]. The situation is even more favourable in panel (b), where the deviations that 
are expected for larger Qr are suppressed by the strong nonlinearity. These numerical results offer a further validation 
of the analytical approximations underlying our the photonic approach. 
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